Our goal is to simulate optical cloaking devices by calculating how light travels through defined boundaries. Optical cloaking refers to the act of making something invisible in some directions by preventing the scattering of light as it hits the boundary. We use boundary integral equation methods to compute the solution in layered boundaries.
In this notebook, we implement in Python an approximation through the use of two main functions:
$\Delta u + k^2 u = 0 \hspace{70mm}$ (Hemholtz Wave Equation)
$u(x) = \int_{B}\frac{\partial \Phi (x,y)}{\partial \nu (y)}\mu(y) - ik\int_{B} \Phi (x,y) \mu (y) dy \hspace{22mm} $(Kress's Boundary Integral Equation)
$\Phi$ represents the fundamental solution to the Hemholtz Wave Equation:
$\Phi (x,y) = \frac{i}{4} H_{0}^{(1)}(k\left|x - y\right |)$.
For transmission, we need to extend this method to include the exterior and the interior of the boundary together. Instead of a single boundary integral equation, we use two to define the solution inside and outside.
$u_e(x) = f(x) + \int_{B}\frac{\partial \Phi ^{e}}{\partial \nu(y)}(x,y) u_e(y) - \int_{B}\Phi ^{e}(x,y)\frac{\partial u_e (y)}{\partial \nu (y)} d\sigma _y \quad x$ $\epsilon$ $E$
$u_i(x) = -\int_{B}\frac{\partial \Phi ^{i}}{\partial \nu (y)}(x,y)u_i(y) + \int_{B}\Phi ^{i} (x,y) \frac{\partial u_i(y)}{\partial \nu(y)}(y)d\sigma _y \quad x$ $\epsilon$ $I$
B refers to the domain of the boundary. E and I are the domains outside and inside of the boundary, respectively. Using the transmission conditions, $u_i = u_e$ and $\frac{1}{k_i} \frac{\partial u_i}{\partial \nu} = \frac{1}{k_e} \frac{\partial u_e}{\partial \nu}$, this system of equations can be rewritten in terms of $u_i$ and $\frac{\partial u_i}{\partial \nu}$.
$u_e(x) = f(x) + \int_{B}\frac{\partial \Phi ^{e}}{\partial \nu(y)}(x,y) u_i(y) - \frac{k_e}{k_i}\int_{B}\Phi ^{e}(x,y)\frac{\partial u_i (y)}{\partial \nu (y)} d\sigma _y \quad x$ $\epsilon$ $E$
$u_i(x) = -\int_{B}\frac{\partial \Phi ^{i}}{\partial \nu (y)}(x,y)u_i(y) + \int_{B}\Phi ^{i} (x,y) \frac{\partial u_i(y)}{\partial \nu(y)}(y)d\sigma _y \quad x$ $\epsilon$ $I$
These equations can be rewritten and turned into a matrix equation. By solving this matrix system for $u_i$ and $\frac{\partial u_i}{\partial \nu}$ on the boundary, those values can be used to approximate via the PTR inside and outside the boundary.
$$\begin{bmatrix} \frac{I}{2} - D^{e} & \frac{k_e}{k_i}S^{e} \\ \frac{I}{2} + D^{i} & -S^{i} \end{bmatrix} \begin{bmatrix}u_i \\ \frac{\partial u_i}{\partial \nu}\end{bmatrix} = \begin{bmatrix} f \\ 0 \end{bmatrix}$$
class Boundary:
'''
Defines the class Boundary, that contains geometric features of a boundary
Attributes
==========
param: parametrization of a boundary
veloc: derivative of the parametrization
accel: second derivative of the parametrization
normal: normal vector to the boundary
tangent: tangent vectors to the boundary
curvature: signed curvature of the boundary
jacobian: Jacobian of the boundary
'''
def __init__(self, b):
self.y = b
@property
def yp(self):
return [sp.diff(p, t) for p in self.y]
@property
def ypp(self):
return [sp.diff(p, t) for p in self.yp]
@property
def J(self):
v2 = [i**2 for i in self.yp]
sv2 = sum(v2)
print(sv2) # This is correct ...
return sp.sqrt(sv2)
@property
def Ï„(self):
return [p/self.J for p in self.yp]
@property
def ν(self):
if len(self.y) == 2:
tmp = (self.yp[1]/self.J, -self.yp[0]/self.J )
return tmp
else:
raise ValueError('Need to define the normal vector for higher dimensions')
@property
def κ(self): #curvature kappa
if len(self.y) == 2:
tmp = (self.yp[0]*self.ypp[1] - self.yp[1]*self.ypp[0])/self.J**3
return tmp
else:
raise ValueError('Need to define the mean curvature for higher dimensions')
def find_boundary_data(f, xbdy, ybdy, ν, J, κ):
'''
Computes u and du/dv on the boundary by solving the BIE system
Attributes
======
kernel of the double layer, L
kernel of the single layer, M
solution on the boundary, u
normal derivative on the boundary, dvu
'''
E_C = 0.5772156649015328606065120900824; # Euler's constant
N = M / 2;
m = np.arange(1, N);
# Create array for ifft...
a = [0]
a.extend(1/m)
a.append(1/N)
a.extend((1/m)[::-1])
Rj = -2 * np.pi * np.fft.ifft(a);
R = np.real(toeplitz(Rj, Rj)); #check coefficients Matlab Toeplitz, order correct ?
# Prepare arrays for solving integral equation of the boundary (Kress pt. 2)
L_e = np.zeros((M, M), dtype=complex);
M_e = np.zeros((M, M), dtype=complex);
L_i = np.zeros((M, M), dtype=complex);
M_i = np.zeros((M, M), dtype=complex);
L1_e = np.zeros((M, M), dtype=complex);
L2_e = np.zeros((M, M), dtype=complex);
M1_e = np.zeros((M, M), dtype=complex);
M2_e = np.zeros((M, M), dtype=complex);
L1_i = np.zeros((M, M), dtype=complex);
L2_i = np.zeros((M, M), dtype=complex);
M1_i = np.zeros((M, M), dtype=complex);
M2_i = np.zeros((M, M), dtype=complex);
A = np.zeros((2*M, 2*M), dtype=complex);
F = np.zeros(2*M, dtype=complex);
for m in range(0, M):
for n in range(0, M):
rdiff = np.asarray([xbdy[m] - xbdy[n], ybdy[m] - ybdy[n]]);
distance = sqrt(np.power(rdiff[:][0], 2) + np.power(rdiff[:][1], 2));
M1_e[m][n] = -0.5 / np.pi * J[n] * besselj(0, k_e * distance);
M1_i[m][n] = -0.5 / np.pi * J[n] * besselj(0, k_i * distance);
if m == n:
#Implementration of the Kress Quadrature
L1_e[m][n] = 0;
L2_e[m][n] = 0.5 / np.pi * κ[n] * J[n];
M2_e[m][n] = J[n] * ( 0.5 * 1j - E_C / np.pi - 0.5 / np.pi * np.log( 0.25 * np.power(k_e, 2) * np.power(J[n], 2)));
L1_i[m][n] = 0;
L2_i[m][n] = 0.5 / np.pi * κ[n] * J[n];
M2_i[m][n] = J[n] * ( 0.5 * 1j - E_C / np.pi - 0.5 / np.pi * np.log( 0.25 * np.power(k_i, 2) * np.power(J[n], 2)));
else:
#Compute cosθ and logterm for kernel calculation
a = np.asmatrix([ν[0][n], ν[1][n]]);
b = np.asmatrix(rdiff/distance);
cosθ = np.asscalar(np.matmul(a, b.transpose()));
logterm = log(float(4 * np.power(sin(0.5 * (m - n) * dt), 2))); #why float ?
#Compute kernels for the BIE system
L1_e[m][n] = 0.5 * k_e / np.pi * J[n] * cosθ * besselj(1, k_e * distance);
L2_e[m][n] = 0.5 * 1j * k_e * J[n] * cosθ * hankel1(1, k_e * distance) - L1_e[m][n] * logterm;
M2_e[m][n] = 0.5 * 1j * J[n] * hankel1(0, k_e * distance) - M1_e[m][n] * logterm;
L1_i[m][n] = 0.5 * k_i / np.pi * J[n] * cosθ * besselj(1, k_i * distance);
L2_i[m][n] = 0.5 * 1j * k_i * J[n] * cosθ * hankel1(1, k_i * distance) - L1_i[m][n] * logterm;
M2_i[m][n] = 0.5 * 1j * J[n] * hankel1(0, k_i * distance) - M1_i[m][n] * logterm;
L_e[m][n] = 0.5 * R[m][n] * (L1_e[m][n]) + 0.5 * (np.pi / N) * (L2_e[m][n]);
M_e[m][n] = (0.5 * R[m][n] * (M1_e[m][n]) + 0.5 * (np.pi / N) * (M2_e[m][n]));
L_i[m][n] = 0.5 * R[m][n] * (L1_i[m][n]) + 0.5 * (np.pi / N) * (L2_i[m][n]);
M_i[m][n] = 0.5 * R[m][n] * (M1_i[m][n]) + 0.5 * (np.pi / N) * (M2_i[m][n]);
#Matrix of combined representation formulas
A11 = 0.5 * np.identity(M) - L_e;
A12 = k_ratio * M_e;
A21 = 0.5 * np.identity(M) + L_i;
A22 = -M_i;
#Combine A's into one matrix, A
A[0:M, 0:M] = A11; A[0:M, M: 2*M] = A12;
A[M: 2*M, 0:M] = A21; A[M: 2*M, M: 2*M] = A22;
#Make matrix for f and 0's, call it F
F[0:M] = f;
#Solve for u and dvu from AU = F
U = np.linalg.solve(A,F);
u = U[0:M];
dvu = U[M: 2*M];
return u, dvu
def make_solution_grid(ngrid, f, xbdy, ybdy, ν, J, κ):
'''
Calculate the solution for interior and exterior domains using u and dvu
Attributes
==========
domain and range values for plot of the solution grid, x and y
approximation of the solution using PTR, v
exact solution for known case of ke = ki, exact
'''
# Find the absolute maximum curvature |κ|max
κ_max = np.amax(np.abs(κ));
dn = (1/κ_max)/ngrid;
rgrid = np.arange(0+dn, 1/ κ_max, dn);
# Allocate memory for the solution
x_e = np.zeros((M + 1, ngrid - 1));
y_e = np.zeros((M + 1, ngrid - 1));
v_e = np.zeros((M + 1, ngrid - 1), dtype=complex); # Kernel and sum will include complex numbers
exact_e = np.zeros((M + 1, ngrid - 1), dtype=complex); # Kernel and sum will include complex numbers
x_i = np.zeros((M + 1, ngrid - 1));
y_i = np.zeros((M + 1, ngrid - 1));
v_i = np.zeros((M + 1, ngrid - 1), dtype=complex); # Kernel and sum will include complex numbers
exact_i = np.zeros((M + 1, ngrid - 1), dtype=complex); # Kernel and sum will include complex numbers
for m in range(0, M):
for n in range(0, ngrid-1):
x_e[m][n] = xbdy[m] + rgrid[n]*ν[0][m];
y_e[m][n] = ybdy[m] + rgrid[n]*ν[1][m];
xe_diff = x_e[m][n] - xbdy;
ye_diff = y_e[m][n] - ybdy;
distance_e = sqrt(np.power(xe_diff, 2) + np.power(ye_diff, 2)); #(r in kress)
cosθ_e = (ν[0][:]*xe_diff + ν[1][:]*ye_diff ) / distance_e;
kernelD_e = 0.25 * 1j * (k_e * cosθ_e * hankel1(1, k_e*distance_e));
kernelS_e = 0.25 * 1j * hankel1(0, k_e*distance_e); # kress (2.2+)
fx = np.exp(1j * k_e * (np.cos(alpha) * x_e[m][n] + np.sin(alpha) * y_e[m][n]));
v_e[m][n] = fx + sum(kernelD_e * J * u)* dt - k_ratio * sum(kernelS_e * J * dvu) * dt; #From Representation Formula (1)
exact_e[m][n] = fx;
# Do the same as above but for the interior...
x_i[m][n] = xbdy[m] - rgrid[n]*ν[0][m];
y_i[m][n] = ybdy[m] - rgrid[n]*ν[1][m];
xi_diff = x_i[m][n] - xbdy;
yi_diff = y_i[m][n] - ybdy;
distance_i = sqrt(np.power(xi_diff, 2) + np.power(yi_diff, 2)); #(r in kress)
cosθ_i = (ν[0][:]*xi_diff + ν[1][:]*yi_diff) / distance_i;
kernelD_i = 0.25 * 1j * (k_i * cosθ_i * hankel1(1, k_i*distance_i));
kernelS_i = 0.25 * 1j * hankel1(0, k_i*distance_i);
v_i[m][n] = -sum(kernelD_i * J * u) * dt + sum(kernelS_i * J * dvu) * dt; #From Representation Formula (2)
exact_i[m][n] = np.exp(1j * k_i * (np.cos(alpha) * x_i[m][n] + np.sin(alpha) * y_i[m][n]));
#Impose Periodicity
x_e[M][:] = x_e[0][:];
y_e[M][:] = y_e[0][:];
v_e[M][:] = v_e[0][:];
exact_e[M][:] = exact_e[0][:];
x_i[M][:] = x_i[0][:];
y_i[M][:] = y_i[0][:];
v_i[M][:] = v_i[0][:];
exact_i[M][:] = exact_i[0][:];
#Combine results for plot displays
x = np.zeros((2 * (M + 1), ngrid - 1));
y = np.zeros((2 * (M + 1), ngrid - 1));
v = np.zeros((2 * (M + 1), ngrid - 1), dtype=complex); # Kernel and sum will include complex numbers
exact = np.zeros((2 * (M + 1), ngrid - 1), dtype=complex); # Kernel and sum will include complex numbers
x[0:M+1][:] = x_i; x[M+1:2 * (M+1)] = x_e;
y[0:M+1][:] = y_i; y[M+1:2 * (M+1)] = y_e;
v[0:M+1][:] = v_i; v[M+1:2 * (M+1)] = v_e;
exact[0:M+1][:] = exact_i; exact[M+1:2 * (M+1)] = exact_e;
return x, y, v, exact
In order to validate our method, we will firstly display the results for a case such that the solution is known. This solution is such that:
$$f(x) = e^{i k_e (cos(\alpha)x + sin(\alpha)y)}$$
$$k_e = k_i$$
The values of $\alpha$, $k_e$, and $k_i$ are chosen below when we set the paramaters for our calculations. Our solution is given by the value $v$ and is compared to the known solution, which we call $exact$. This error is displayed on a log plot which follows the surface plots of the approximate and exact solutions.
Our method is then use to examine our results for different scenarios. For example, we apply our method for the case when $k_i = -0.5 k_e$. We also test our method for different shapes of the boundary, such as a star-shaped boundary.
# Specify the shape of the boundaries
t = sp.Symbol('t')
Ellipse = (2*cos(t), 3*sin(t))
Circle = (cos(t), sin(t))
Star = ((1+0.3*cos(5*t))*cos(t), (1+0.3*cos(5*t))*sin(t))
Kite = ((cos(t) + 0.65*cos(2*t) -0.65, 1.5*sin(t)))
E = Boundary(Ellipse)
K = Boundary(Kite)
C = Boundary(Circle)
S = Boundary(Star)
# Convert the sympy into numpy quantities
bdy_E = sp.lambdify(t, E.y);
ν_E = sp.lambdify(t, E.ν);
J_E = sp.lambdify(t, E.J);
κ_E = sp.lambdify(t, E.κ);
bdy_K = sp.lambdify(t, K.y);
ν_K = sp.lambdify(t, K.ν);
J_K = sp.lambdify(t, K.J);
κ_K = sp.lambdify(t, K.κ);
bdy_C = sp.lambdify(t, C.y);
ν_C = sp.lambdify(t, C.ν);
J_C = sp.lambdify(t, C.J);
κ_C = sp.lambdify(t, C.κ);
bdy_S = sp.lambdify(t, S.y);
ν_S = sp.lambdify(t, S.ν);
J_S = sp.lambdify(t, S.J);
κ_S = sp.lambdify(t, S.κ);
4*sin(t)**2 + 9*cos(t)**2 4*sin(t)**2 + 9*cos(t)**2 4*sin(t)**2 + 9*cos(t)**2 4*sin(t)**2 + 9*cos(t)**2 (-sin(t) - 1.3*sin(2*t))**2 + 2.25*cos(t)**2 (-sin(t) - 1.3*sin(2*t))**2 + 2.25*cos(t)**2 (-sin(t) - 1.3*sin(2*t))**2 + 2.25*cos(t)**2 (-sin(t) - 1.3*sin(2*t))**2 + 2.25*cos(t)**2 sin(t)**2 + cos(t)**2 sin(t)**2 + cos(t)**2 sin(t)**2 + cos(t)**2 sin(t)**2 + cos(t)**2 (-(0.3*cos(5*t) + 1)*sin(t) - 1.5*sin(5*t)*cos(t))**2 + ((0.3*cos(5*t) + 1)*cos(t) - 1.5*sin(t)*sin(5*t))**2 (-(0.3*cos(5*t) + 1)*sin(t) - 1.5*sin(5*t)*cos(t))**2 + ((0.3*cos(5*t) + 1)*cos(t) - 1.5*sin(t)*sin(5*t))**2 (-(0.3*cos(5*t) + 1)*sin(t) - 1.5*sin(5*t)*cos(t))**2 + ((0.3*cos(5*t) + 1)*cos(t) - 1.5*sin(t)*sin(5*t))**2 (-(0.3*cos(5*t) + 1)*sin(t) - 1.5*sin(5*t)*cos(t))**2 + ((0.3*cos(5*t) + 1)*cos(t) - 1.5*sin(t)*sin(5*t))**2
# Evaluate boundary data
def def_f(bdy):
xbdy = bdy[0][:];
ybdy = bdy[1][:];
# Plane Wave Equation
f = np.exp(1j * k_e * (np.cos(alpha) * xbdy + np.sin(alpha) * ybdy));
return xbdy, ybdy, f
# Evaluation of the results for given parameters
# USER INPUTS: Modify M, k_e, k_i, and ngrid to modify results
M = 300
k_e = 3
k_i = 3
k_ratio = k_e/k_i;
ngrid = 49;
alpha = np.pi / 5;
# Values for parameterization of the boundaryh
dt = 2.0 * np.pi/(M);
θ = np.arange(0, 2* np.pi, dt);
# Using bdy, ν, J, and κ function names for ellipse
xbdy, ybdy, f = def_f(bdy_E(θ));
u, dvu = find_boundary_data(f, xbdy, ybdy, ν_E(θ), J_E(θ), κ_E(θ));
df = [1j * k_e * np.cos(alpha) * np.exp(1j * k_e * (np.cos(alpha) * xbdy + np.sin(alpha) * ybdy)), 1j * k_e * np.sin(alpha) * np.exp(1j * k_e * (np.cos(alpha) * xbdy + np.sin(alpha) * ybdy))]
dvf = ν_E(θ)[0] * df[0] + ν_E(θ)[1] * df[1]
fig, ax = plt.subplots(1,1)
fig.suptitle('Boundary data comparision');
ax.plot(np.real(u), 'b', label='Re(f)')
ax.plot(np.real(f), 'r--', label='Re(u)')
ax.plot(np.imag(u), 'k', label='Im(f)')
ax.plot(np.imag(f), 'm--', label='Im(u)')
ax.legend(loc='upper right');
fig, ax = plt.subplots(1,1)
ax.plot(np.real(dvu), 'b', label='Re(dvf)')
ax.plot(np.real(dvf), 'r--', label='Re(dvu)')
ax.plot(np.imag(dvu), 'k', label='Im(dvf)')
ax.plot(np.imag(dvf), 'm--', label='Im(dvu)')
ax.legend(loc='upper right');
#TO DO check dvu as well
x, y, v, exact = make_solution_grid(ngrid, f, xbdy, ybdy, ν_E(θ), J_E(θ), κ_E(θ));
fig = plt.figure(figsize=(15,12))
fig.suptitle("Ellipse Boundary for k_e = k_i");
ax = fig.add_subplot(2 ,2 ,1, projection='3d')
ax.view_init(90,90)
real1 = ax.plot_surface(x, y, np.real(v), rstride=1, cstride=1, cmap=plt.get_cmap('viridis'), linewidth=0, antialiased=False)
plt.title("Real Part of the PTR Solution");
fig.colorbar(real1, shrink=0.5);
plt.axis('off')
plt.grid(b=None)
ax = fig.add_subplot(2, 2, 2, projection='3d')
ax.view_init(90,90)
imag1 = ax.plot_surface(x, y, np.imag(v), rstride=1, cstride=1, cmap=plt.get_cmap('viridis'), linewidth=0, antialiased=False)
plt.title("Imaginary Part of the PTR Solution");
fig.colorbar(imag1, shrink=0.5);
plt.axis('off')
plt.grid(b=None)
ax = fig.add_subplot(2 ,2 ,3, projection='3d')
ax.view_init(90,90)
real1 = ax.plot_surface(x, y, np.real(exact), rstride=1, cstride=1, cmap=plt.get_cmap('viridis'), linewidth=0, antialiased=False)
plt.title("Real Part of the Exact Solution");
fig.colorbar(real1, shrink=0.5);
plt.axis('off')
plt.grid(b=None)
ax = fig.add_subplot(2, 2, 4, projection='3d')
ax.view_init(90,90)
imag1 = ax.plot_surface(x, y, np.imag(exact), rstride=1, cstride=1, cmap=plt.get_cmap('viridis'), linewidth=0, antialiased=False)
plt.title("Imaginary Part of the Exact Solution");
fig.colorbar(imag1, shrink=0.5);
plt.axis('off')
plt.grid(b=None)
data= [Surface(x=x, y=y, z=log(np.absolute(v - exact)))]
layout=dict(xaxis=dict(zeroline=False),
yaxis=dict(zeroline=False),
width=600,
height=600,
margin=dict(l=50, r=50, b=100, t=100),
title='Log Error Plot', hovermode='closest'
)
figure1=dict(data=data, layout=layout)
iplot(figure1)
# TO DO: isseu with plotly at (x,y) = (pi,0) ? To investigate
M = 300
k_e = 2
k_i = 4
k_ratio = k_e/k_i;
ngrid = 49;
# Values for parameterization of the boundaryh
dt = 2.0 * np.pi/(M);
θ = np.arange(0, 2* np.pi, dt);
# Using bdy, ν, J, and κ function names for ellipse
xbdy2, ybdy2, f2 = def_f(bdy_E(θ));
u, dvu = find_boundary_data(f2, xbdy2, ybdy2, ν_E(θ), J_E(θ), κ_E(θ));
x2, y2, v2, exact2 = make_solution_grid(ngrid, f2, xbdy2, ybdy2, ν_E(θ), J_E(θ), κ_E(θ));
fig = plt.figure(figsize=(15,7))
fig.suptitle("Ellipse Boundary for k_e = 0.5k_i");
ax = fig.add_subplot(1 ,2 ,1, projection='3d')
ax.view_init(90, 90)
real1 = ax.plot_surface(x2, y2, np.real(v2), rstride=1, cstride=1, cmap=plt.get_cmap('viridis'), linewidth=0, antialiased=False)
plt.title("Real Part of the PTR Solution");
fig.colorbar(real1, shrink=0.5);
plt.axis('off')
plt.grid(b=None)
ax = fig.add_subplot(1, 2, 2, projection='3d')
ax.view_init(90,90)
imag1 = ax.plot_surface(x2, y2, np.imag(v2), rstride=1, cstride=1, cmap=plt.get_cmap('viridis'), linewidth=0, antialiased=False)
plt.title("Imaginary Part of the PTR Solution");
fig.colorbar(imag1, shrink=0.5);
plt.axis('off')
plt.grid(b=None)
data= [Surface(x=x2, y=y2, z=np.imag(v2))]
layout=dict(xaxis=dict(zeroline=False),
yaxis=dict(zeroline=False),
width=600,
height=600,
margin=dict(l=50, r=50, b=100, t=100),
title='Movable Plot of the Real Part', hovermode='closest'
)
figure2=dict(data=data, layout=layout)
iplot(figure2)
M = 300
k_e = 2
k_i = -4
k_ratio = k_e/k_i;
ngrid = 49;
# Values for parameterization of the boundaryh
dt = 2.0 * np.pi/(M);
θ = np.arange(0, 2* np.pi, dt);
# Using bdy, ν, J, and κ function names for ellipse
xbdy3, ybdy3, f3 = def_f(bdy_E(θ));
u, dvu = find_boundary_data(f3, xbdy3, ybdy3, ν_E(θ), J_E(θ), κ_E(θ));
x3, y3, v3, exact3 = make_solution_grid(ngrid, f3, xbdy3, ybdy3, ν_E(θ), J_E(θ), κ_E(θ));
fig = plt.figure(figsize=(15,7))
fig.suptitle("Ellipse Boundary for k_e = -0.5k_i");
ax = fig.add_subplot(1 ,2 ,1, projection='3d')
ax.view_init(90,90)
real1 = ax.plot_surface(x3, y3, np.real(v3), rstride=1, cstride=1, cmap=plt.get_cmap('viridis'), linewidth=0, antialiased=False)
plt.title("Real Part of the PTR Solution");
fig.colorbar(real1, shrink=0.5);
plt.axis('off')
plt.grid(b=None)
ax = fig.add_subplot(1, 2, 2, projection='3d')
ax.view_init(90,90)
imag1 = ax.plot_surface(x3, y3, np.imag(v3), rstride=1, cstride=1, cmap=plt.get_cmap('viridis'), linewidth=0, antialiased=False)
plt.title("Imaginary Part of the PTR Solution");
fig.colorbar(imag1, shrink=0.5);
plt.axis('off')
plt.grid(b=None)
data= [Surface(x=x3, y=y3, z=np.imag(v3))]
layout=dict(xaxis=dict(zeroline=False),
yaxis=dict(zeroline=False),
width=600,
height=600,
margin=dict(l=50, r=50, b=100, t=100),
title='Movable Plot of the Real Part', hovermode='closest'
)
figure3=dict(data=data, layout=layout)
iplot(figure3)
#TO DO: create widgets to do displya for all figures
M = 300
k_e = 4
k_i = 4
k_ratio = k_e/k_i;
ngrid = 49;
# Values for parameterization of the boundaryh
dt = 2.0 * np.pi/(M);
θ = np.arange(0, 2* np.pi, dt);
# Using bdy, ν, J, and κ function names for ellipse
xbdy4, ybdy4, f4 = def_f(bdy_S(θ));
u, dvu = find_boundary_data(f4, xbdy4, ybdy4, ν_S(θ), J_S(θ), κ_S(θ));
x4, y4, v4, exact4 = make_solution_grid(ngrid, f4, xbdy4, ybdy4, ν_S(θ), J_S(θ), κ_S(θ));
#TO DO: create the exact solution not in the make_solution_grid
fig = plt.figure(figsize=(15,12))
fig.suptitle("Ellipse Boundary for k_e = k_i");
ax = fig.add_subplot(2 ,2 ,1, projection='3d')
ax.view_init(90,90)
real1 = ax.plot_surface(x4, y4, np.real(v4), rstride=1, cstride=1, cmap=plt.get_cmap('viridis'), linewidth=0, antialiased=False)
plt.title("Real Part of the PTR Solution");
fig.colorbar(real1, shrink=0.5);
plt.axis('off')
plt.grid(b=None)
ax = fig.add_subplot(2, 2, 2, projection='3d')
ax.view_init(90,90)
imag1 = ax.plot_surface(x4, y4, np.imag(v4), rstride=1, cstride=1, cmap=plt.get_cmap('viridis'), linewidth=0, antialiased=False)
plt.title("Imaginary Part of the PTR Solution");
fig.colorbar(imag1, shrink=0.5);
plt.axis('off')
plt.grid(b=None)
ax = fig.add_subplot(2 ,2 ,3, projection='3d')
ax.view_init(90,90)
real1 = ax.plot_surface(x4, y4, np.real(exact4), rstride=1, cstride=1, cmap=plt.get_cmap('viridis'), linewidth=0, antialiased=False)
plt.title("Real Part of the Exact Solution");
fig.colorbar(real1, shrink=0.5);
plt.axis('off')
plt.grid(b=None)
ax = fig.add_subplot(2, 2, 4, projection='3d')
ax.view_init(90,90)
imag1 = ax.plot_surface(x4, y4, np.imag(exact4), rstride=1, cstride=1, cmap=plt.get_cmap('viridis'), linewidth=0, antialiased=False)
plt.title("Imaginary Part of the Exact Solution");
fig.colorbar(imag1, shrink=0.5);
plt.axis('off')
plt.grid(b=None)
data= [Surface(x=x4, y=y4, z=log(np.absolute((v4) - (exact4))))]
layout=dict(xaxis=dict(zeroline=False),
yaxis=dict(zeroline=False),
width=700,
height=700,
margin=dict(l=50, r=50, b=100, t=100),
title='Log Error Plot', hovermode='closest'
)
figure1=dict(data=data, layout=layout)
iplot(figure1)
# TO DO: isseu with plotly at (x,y) = (pi,0) ? To investigate